
function [frame] = SectionPlotsubroutine(setasfile,modelvar)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS FUNCTION GENERATES THE FIGURE OF A SETAS MODEL VARIABLE AND RETURNS 
% THE FRAME FOR THE ANIMATION PRODUCTION.
% 
%
% INPUT  : - setasfile : a setas netcdf file name
%          - modelvar  : a SETAS variable. Available options are 'temp or 
%                       'salt' 
% OUTPUT :  frame : image frame  
%
% Author: Benedicte Pasquer, IMOS/eMII (http://imos.org.au/)
% email: benedicte.pasquer@utas.edu.au
% May 2013   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%section path coordinates
section = load('plum_Derwentsect1.sect');
xsectc  = section(:,1);
ysectc  = section(:,2);
%SERIES OF WAYPOINTS ALONG THE SECTION PATH
pointsxy = load('plum_Derwentsect.xy');
npoints  = length(xsectc) ;

%EXTRACT DATA FROM SETAS NETCDF FILE
ncid = netcdf.open(setasfile,'NC_NOWRITE');
 
varid = netcdf.inqVarID(ncid,'latitude');
lat = netcdf.getVar(ncid,varid); 
varid = netcdf.inqVarID(ncid,'longitude');
lon = netcdf.getVar(ncid,varid); 
varid = netcdf.inqVarID(ncid,'zc') ;
z_centre = netcdf.getVar(ncid,varid);
varid = netcdf.inqVarID(ncid,'botz') ;
botz = netcdf.getVar(ncid,varid);
%GET SPECIFIED VARIABLE DATA
switch modelvar
    case 'temp'
        Var4title ='Temperature [degC]'; 
        varid = netcdf.inqVarID(ncid,'temp');
        
    case 'salt'
        Var4title ='Salinity [ppt]';
        varid = netcdf.inqVarID(ncid,'salt');

end

var = netcdf.getVar(ncid,varid,'double');
var = mean(var,4);
var = squeeze(var);%REMOVE THE SINGLETON DIMENSION    
var = permute(var,[2 1 3]);
netcdf.close(ncid)

%DETERMINE TOP K LEVEL IN FILE & READ IN DEPTH LEVEL.

nk = length(z_centre) ;
z_grid(2:nk) = 0.5*(z_centre(1:nk-1) + z_centre(2:nk)) ;
z_grid(nk+1) = z_grid(nk) + 2*(z_centre(nk)- z_grid(nk)) ;
z_grid(1) = z_grid(2) - 2*(z_grid(2) - z_centre(1)) ;


%HORIZONTAL GRID GEOMETRY.
%DETERMINE GRID CELL COORDS FOR THE SECTION CELL CENTRES
 
if (~isempty(lon)),

    lon = lon' ; lat = lat' ;
    x_centre = lon ; y_centre = lat ;
    
end
%DETERMINE CELL CENTRE LOCATION
[m,n] = size(x_centre) ;
x0 = 2*x_centre(1,:) - x_centre(2,:) ;
y0 = 2*y_centre(1,:) - y_centre(2,:) ;
xf = 2*x_centre(m,:) - x_centre(m-1,:) ;
yf = 2*y_centre(m,:) - y_centre(m-1,:) ;
xs = ([x0 ; x_centre] + [x_centre ; xf])/2 ;
ys = ([y0 ; y_centre] + [y_centre ; yf])/2 ;
x0 = 2*xs(:,1) - xs(:,2) ;
y0 = 2*ys(:,1) - ys(:,2) ;
xf = 2*xs(:,n) - xs(:,n-1) ;
yf = 2*ys(:,n) - ys(:,n-1) ;
x_grid = ([x0 xs] + [xs xf])/2 ;
y_grid = ([y0 ys] + [ys yf])/2 ;

% DETERMINE HORIZONTAL GRID CELL COORDS. FOR THE SECTION CELL CENTRES

[isect,jsect] = xy2ij_curv (xsectc,ysectc,x_grid,y_grid,2) ;

% DETERMINE DEPTH VALUES ALONG SECTION.
botz = botz';
botz= -botz ;
for i=1:length(section)
    
     var_sect(1:nk,i) = squeeze(var(jsect(i),isect(i),:));
     bot_sect(i) = botz(jsect(i),isect(i));

end

var_sect(var_sect>1000) = NaN; 
bot_sect(bot_sect==99) = NaN;   
 
% PREPARE THE GRID FOR PCOLOR. REPLACE GRID CELL OF BOTTOM BY BOTTOM  VALUE
% ADD A RAW OF DUMMY DATA TO GET THE PCOLOR TO WORK PROPERLY 
var_sect(23,:) = NaN;

% COMPUTES DISTANCE BETWEEN SECTION POINTS
for  j = 2 : length(section) 
    dd(j) = distance(section(j-1,2),section(j-1,1),section(j,2),section(j,1));
    dd(j) = deg2km(dd(j));
    
end
sitedist = cumsum(dd(1:end));

[X,Y] = meshgrid(sitedist,z_grid);
for i = 1:npoints
    ind = find(isnan(var_sect(:,i)),2,'last');
    Y(ind(1)+1,i) = bot_sect(i);
end

% ADD LOCATION OF SECTION SITES ON THE FIGURE
% COMPUTES DISTANCE BEWEEN POINTS
for  jp = 2 : length(pointsxy) 
    ddp(jp) = distance(pointsxy(jp-1,2),pointsxy(jp-1,1),pointsxy(jp,2),pointsxy(jp,1));
    ddp(jp) = deg2km(ddp(jp));
    
end
sitedistp = cumsum(ddp(1:end));

% FIGURE
hfr = figure('Color',[.9 .98 .98]);
axes1 = axes('Parent',hfr,'XAxisLocation','top',...
    'Layer','top',...
    'CLim',[min(min(var_sect)) max(max(var_sect))],...
    'Color', [.9 .98 .98],'YColor', [0 0 .7], 'Xcolor', [0 0 .7],...
    'Fontname','arial','Fontsize',12,'Fontweight','bold','box','on');

box(axes1,'on');
hold(axes1,'all');
% PCOLOR
pcolor(X,Y,var_sect,'Parent',axes1)
shading flat
axis tight

h = line(sitedistp,zeros(1,length(sitedistp)));
set(h,'Marker','.','Markersize',8,'Color','k','Parent',axes1);
line(sitedistp(1),0,'Marker','*','Markersize',8,'Color','r','Parent',axes1);
yLim([-60 0])
ylabel('depth [m]','Fontsize',10)
xlabel('Distance from start [km]','Fontsize',10)

dash = regexp(setasfile,'-');
yr = setasfile(dash-4:dash-1);
mth = setasfile(dash+1:dash+2);

title({['SETAS monthly mean ',Var4title,  ' section ',mth,'/',yr];'along the Derwent estuary'},...
        'Fontname','arial','Fontsize',12,'Fontweight','bold', ...
        'color',[0 0 .7]);
    
% CREATE COLORBAR
if strcmp(Var4title(1:4),'Temp') == 1
    caxis([5 20])
    cbarytick =[5 8 11 14 17 20];
else %Salinity
    caxis([0 35])
     cbarytick =[0 5 10 15 20 25 30 35];
end
colorbar('peer',axes1,'NextPlot','replacechildren',...
    'Ytick',cbarytick,...
    'Fontname','arial','Fontsize',10,'Fontweight','bold', ...
        'Color',[0 0 .7]);

axes2 = axes( 'Position', [0.16 0.185 0.26 0.36]);

% DISPLAY SECTION PATH IN FIGURE INSET 
m_proj('mercator','long',[147.15 147.5],'lat',[-43.1 -42.7])

m_line(xsectc,ysectc,'Linewidth',1,'Color','r','Parent',axes2)
m_line(pointsxy(:,1),pointsxy(:,2),'Marker','.',...
    'Markersize',8,...
    'Color', 'k',...
    'LineStyle','none',...
    'Parent',axes2);
m_line(pointsxy(1,1),pointsxy(1,2),'Marker','*',...
    'Markersize',6,...
    'Color', 'r',...
    'LineStyle','none',...
    'Parent',axes2);

m_usercoast('Derwent.mat','patch',[.9 .9 .9],'edgecolor','none');
m_grid('xtick',[147 147.16667 147.3333 147.5],...
    'ytick',[-43 -42.8333 -42.9166 -42.75],...
    'linewi',1,'linest',':',...
    'tickdir','out','fontsize',8);
title('Section path','Fontsize',8)

%PREPARE FRAME FOR ANIMATION
drawnow
frame = getframe(hfr);

close(gcf)





